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Abstract 

For low-pressure discharges, when the electron mean free path is larger or comparable with 
the discharge length, the electron dynamics is essentially nonlocal. Moreover, the electron energy 
distribution function (EEDF) deviates considerably from a Maxwellian. Therefore, an accurate 
kinetic description of the low-pressure discharges requires knowledge of the nonlocal conductiv- 
ity operator and calculation of the nonMaxwellian EEDF. The previous treatments made use of 
simplifying assumptions: a uniform density profile and a Maxwellian EEDF. In the present study 
a self-consistent system of equations for the kinetic description of nonlocal, nonuniform, nearly 
collisionless plasmas of low-pressure discharges is derived. It consists of the nonlocal conductiv- 
ity operator and the averaged kinetic equation for calculation of the nonMaxwellian EEDF. The 
importance of accounting for the nonuniform plasma density profile on both the current density 
profile and the EEDF is demonstrated. 
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List of variables: 

Tb is the bounce period, 

T is half of the bounce period, T(e x ) = T^/2 = /*+ dxj \v x (x, e x )\, 
fib is the bounce frequency, = 2n/T h , 
L is the gap width, 

R is half of the gap width, which is used to model a cylinder geometry, R = L/2, 

w is the electron kinetic energy, w = w x + w y + w z = m{y 2 x + v 2 + v 2 )/2, 

(f(x) is the electron potential energy, <p(x) = — ecj)(x), and 

(ft is the electrostatic potential, 

e is the total electron energy, e — w + <p, 

uj is the frequency of the rf electric field, 

v is the electron elastic collision frequency, 

A is the electron elastic mean free path, 

z/£ is the electron inelastic collision frequency for the process number k, 

x + (e) are the left and right turning points [e = <p(x±)), 
r is the time of flight from the left turning point x_(e x ) to x: r(x, e x ) = f*_ dxj \v x (x, e x )\, 
9 is the variable angle for bounce motion, defined as 9(x) = nsgn(v x )r(x), /T{e x ) 

e x ) is the generalized phase of the rf electric field, e x ) = ^ {—iuj + 

v)dxj \v x (x,e x )\. If uj » v, $ ~ ur, 
5 is the width of the skin layer, 

/ is the electron velocity distribution function (EVDF), / = fo(e) + fi, where fo(e) is 
the main part of the electron velocity distribution function averaged over velocity directions 
and over available space for electrons with a given total energy e, which is referred in the 
following as the electron energy probability function (EEPF). Notwithstanding the fact that 
/ is defined in the velocity space, we shall look for f(e) as a function of the energy. The EEPF 
is normalized as n = J fd 3 v = 47rv / 2/m 3//2 J^ x ^ fo(e)Je — <p(x)de, where n is the electron 

density, and the factor 4:ir\/2/m 3 ^ 2 is included in the definition of fo for convenience, fi is 
the rapidly varying, anisotropic part of the EVDF, 

E sc (x) is the space charge stationary electric field, 

E y (x,t) is the rf nonstationary electric field, 

St(f) is the collision integral, 

V r f is the oscillatory electron velocity driven by rf electric field. 

I. INTRODUCTION 

Low pressure radio-frequency discharges are extensively utilized for plasma processing 
and lighting |T|]. Simulation of discharge properties is a common tool for optimization of 
the plasma density profiles and ion and electron fluxes. Recent plasma technology tends to 
decrease the gas pressures down to the millitorr range. For these low pressures it is eas- 
ier to maintain uniform plasmas with well controlled parameters. Due to the large value 
of the electron mean free path (A) the electron current is determined not by the local rf 
electric field (Ohm's law), but rather is a function of the whole profile of the rf electric 
field on distances of order A (anomalous skin effect). Therefore, a rather complicated non- 
local conductivity operator has to be determined for the calculation of the rf electric field 
penetration into the plasma. Moreover, the electron energy distribution function (EEDF) 
is typically nonMaxwellian in these discharges ||. Hence, for accurate calculation of the 
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discharge characteristics at low pressures, the EEDF needs to be computed self-consistently. 
Self-consistency is an important and difficult issue for the kinetic simulations of a plasma. 
The EEDF, nonlocal conductivity and plasma density profiles are all nonlinear and nonlo- 
cally coupled. That is why, the self-consistency aspect of the model is the main concern 
of this study. The so-called "nonlocal" approach relies on the direct semi- analytic solution 
of the Boltzmann equation in the limiting regime where the electron relaxation length is 
much large than the discharge gap, but the electron mean free path is small compared with 
the discharge dimension ||, [Q. The nonlocal approach has been successfully applied to 
the self-consistent kinetic modelling of various low-pressure discharges, where the electron 
mean free path is small: the capacitively coupled plasmas ||; |§, the inductively coupled 
|| , lK| ; the dc discharges |~T| , |L2| ; the afterglow [|13|] , and the surface- wave 
The additional references can be found in reviews: fl5|l , fL6| , and fT 



0- 



plasmas 
discharges 

If gas pressure lowered even farther, the electron mean free path becomes comparable or 
even larger than the discharge dimension and numerous collisionless phenomena dominate 
the discharge characteristics |L8 |. In the present paper the nonlocal approach is general- 
ized for the low-pressure discharges to account for the collisionless heating and transit-time 
(electron inertia) effects in the discharge description. 

Present analysis considers only an inductively coupled plasma. But the approach has 
been designed in the most generalized way, so that derivations can be readily performed for 
other discharges. For example, in Ref. |19[ the capacitive discharge; in Ref. |20| the electron- 



cyclotron-resonance discharge and in Ref. |21[ surface-wave discharge were considered with 
self-consistent account for collisionless heating. 

Most previously reported theoretical studies assume a uniform plasma, in a semi-infinite 
Tl\ or a slab geometry [23|. In this case the analytical treatment considerably simplifies, 
because electron trajectories are straight. In the semi-infinite case, electrons traverse the 
region of the rf electric field (skin layer) and are reflected back into the plasma at the dis- 
charge walls. An acquired velocity kick then dissipates in the plasma on distances of order 
the electron mean free path. If the plasma dimension is small or comparable with A, the 
subsequent kicks are correlated. The resonance between the wave frequency and the bounce 
frequency of the electron motion between walls may result in modification of the nonlocal 
conductivity [p5 and may yield an enhanced electron heating [26| . The anomalous skin 
effect has been studied experimentally in cylindrical |24| and planar discharges [27|. Addi- 
tional references can be found in the reviews of classical and recent works on the anomalous 
skin effect in gas discharge plasmas |28| . 



29| . The theoretical studies in cylindrical geometry 



are much more cumbersome, and has been done for uniform plasma in Refs. p0| , [[31], |32| 



and for a parabolic potential well in Ref. ||33|| . Qualitative results in the cylindrical geom- 
etry are similar to the results in the plane geometry, therefore, in the present study only 
one-dimensional slab geometry is considered. 

For the case of a bounded uniform plasma, the electrostatic potential well is flat in the 
plasma and infinite at the wall (to simulate the existence of sheaths). In this square potential 
well, electrons are reflected back into the plasma only at the discharge walls. In a realistic 
non-uniform plasma, however, the position of the turning points will depend on the electron 
total (kinetic plus potential) energy and the actual shape of the potential well, i.e., low total 
energy electrons bounce back at locations within the plasma and may not reach regions of 
high electric field at all. As a result the current density profiles in a nonuniform plasma 
may considerably differ from the profiles in a uniform plasma. The theory of the anomalous 
skin effect for an arbitrary profile of the electrostatic potential and a Maxwellian EVDF was 
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developed by Meierovich et al. in Refs.[|34]], [0, and |36| for the slab geometry. Although 
rigorous analytical results of non-uniform plasmas have been reported, the detailed self- 
consistent simulations related to such plasmas and comparison with experimental data are 
lacking. Such simulations were completed recently and presented in our separate publications 
|J7J and and will be additionally reported elsewhere. 

The kinetic description of the anomalous skin effect is based on a well known mechanism 
of collisionless power dissipation -the Landau damping |5?J. In the infinite plasma, the 
resonance particles moving with a velocity (v) close to the wave phase velocity, so that 
uj = v ■ k, intensively interact with wave fields. Therefore, the collisionless electron heating 
(and the real part of the surface impedance) depends on the magnitude of a Fourier harmonic 
of the electric field [-E(fc)] and the number of the resonant particles [f(v x = u;//c),x = k//c]. 
That is why, the momentum acquired in the skin layer of width S is maximal if the projection 
of velocity perpendicular to the plasma boundary (x - axis direction) is of order lu/5. If the 
interaction with the skin layer are repeated in a resonance manner the momentum changes 
are mounted up. Therefore, the main contribution to the electron heating and the resistive 
part of the surface impedance comes from these resonant electrons. In a bounded plasma, 
the resonance condition requires the bounce period be equal to one or several rf electric 
field periods: Tf, = 2im/u, where n is an integer number. The maximum interaction occurs 
for n = 1 (see below). For a slab of width L, Tf, = 2L/v x . The maximum electron heating 
occurs if both aforementioned conditions are satisfied simultaneously, which gives u/S = v x 
and 2L/v x = 2tc/uj or L = 5tt |25j. Hence, the optimum conditions for the power transfer to 
the plasma corresponds to the plasma of size comparable with the 3 times of the skin depth. 
Because the bounce frequency depends on the electrostatic potential, accounting for the 
plasma nonuniformity is important for a correct calculation of the efficient power coupling. 

As discussed before, the collisionless heating is determined by the number of resonant 
particles, and, hence, is dependent on the EEDF. The EEDF, in its turn, is controlled by the 
collisionless heating. The only particles, which are in resonance with a wave, are heated by 
the collisionless heating. It means that in the regime of the collisionless dissipation, the form 
of the electron energy distribution function is sensitive to the wave spectrum. Therefore, the 
plateau in the EEDF can be formed in the regions of intensive collisionless heating, if the 



wave phase velocities are confined in some interval [40 . The evidences of a plateau formation 



for the capacitive discharge plasma were obtained in Ref. [fll]]. The cold electrons, which are 
trapped in the discharge center, do not reach periphery plasma regions where an intensive 
rf elective field is located, and as a result, these electrons are not heated by the rf electric 
field. The coupling between the EEDF shape and collisionless heating may result in new 



nonlinear phenomena: an explosive generation of the cold electrons [42]. The experimental 



evidences of the collisionless heating on the EEDF were obtained in Ref. ||27|| , p3| , |44|and 
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In the linear approximation the collisionless dissipation does not depend explicitly on the 
collision frequency. However, as shown in Ref. ||46|| , if the electron elastic collision frequency 
is too small, heating can actually decrease due to nonlinear effects akin of the nonlinear 
Landau damping. 

The present article presents a self-consistent system of equations describing the non-local 
electron kinetics in a 1-D slab (bounded) non-uniform plasma. The system consists of a 
nonlocal conductivity operator, and an averaged over fast electron motions kinetic equation 
for the EEDF. Transit time (non-local) effects on the current density profile and collision- 
less heating are of particular interest. Rigorous derivations for the nonlocal conductivity 
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operator have been performed. The analytic results of Ref.|34| for the Maxwellian EEDF 
were generalized for the nonMaxwellian EEDF. The spectral method was developed to find 
the rf electric field profile. A quasilinear approach was used for calculating the collisionless 
heating. The quasilinear theory developed in Ref . g7| was generalized for an arbitrary value 
of the collision frequency. As a result, the simulations can be done in a wide range of the 
background gas pressures ranging from the collisional case (A << 5) to the fully collisionless 
case (A > L). Self-consistency of the nonlocal conductivity operator and the energy diffusion 
coefficient has been verified: both yield the same expression for the power deposition. The 
robust time- averaging procedure was designed for the kinetic equation in a most general 
way. As a result, the procedure can be readily repeated for other discharges, see for example 
El and |0|. 



II. CALCULATION OF ANISOTROPIC PART OF THE ELECTRON VELOCITY 
DISTRIBUTION FUNCTION f x 



fo + /i, where / (e) 

,2> 



ni\vt + v% + v 2 z )/2 is the kinetic 



In low-pressure discharges, where the energy relaxation length is large compared with the 
plasma width, the main of the electron velocity distribution function (EVDF) is a function 
of the total energy only JFE], [IB| and [IT]]. Therefore, we look for / 
is a function of the total energy e, e = w + <p(x), w = m{y1 + v 7 y , 
energy, (p = —e<ft is the electron electrostatic potential energy, and <ft is the electrostatic 
potential, fi does not contribute to the electron density (the integral f\ over the velocity 
space is equal to zero / fid 3 v = 0), but f\ contributes to the electron current (the integral 
fo over the velocity space weighted with the electron velocity is equal to zero J fovd 3 v = 0). 
Typically the mean electron flow velocity (V r f = J vf\d 3 v/ J fod 3 v) is small compared with 



the thermal velocity (Vr = y2T/m). Therefore, the isotropic part of the EVDF is larger 

than the anisotropic part f 1 ~ (^ r ^/Vr)/o << fo HU; EH' an< ^ 0- 
Vlasov's equation reads: 



dfi ^ dh eE sc (x) df\ eE y (x, t) d(fo + fx) 



Of 



dx 



m 



dv. 



m 



dv 1t 



St{h + f ) 



(2.1) 



where E sc (x) is the space-charge stationary electric field, and E y (x, t) is the rf nonstationary 
electric field, St(f) is the collision integral. In the Eq. (|2.1|) , we used the fact that 



df (e) , eE sc {x) df (e) 



dfo(e) 



dx m dv x dx 



0. 



(2.2) 



because e x is constant along a trajectory. After applying the standard quasilinear theory, 
Eq.( [2.1|) splits into two equations Wh '- a linear equation for f\ 



dt dx 
and a quasilinear equation for f Q 



dfi dfx | eE sc {x) df\ eE y {x , t) df 



ni 



dv„ 



ni 



dv,, 



(2.3) 



eE y (x,t) dfx 



dv,, 



St(fo), 



(2.4) 



where upper bar denotes space-time averaging over the phase space available for the electron 
with the total energy e gt], fi9|,[[50[ ■ 
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The rf electric field E y (x,t) = E y o(x) exp(— iuot) and the anisotropic part of the EVDF 
fi — fio ex p( — iwt) are harmonic functions, where u is the discharge frequency. In what 
follows the subscript is omitted. Eq.([2.3|) becomes 



+ ev y E y (x 



dfo 
de 



(2.5) 



In transformation from Eq. ( |2.3|) to Eq. (|2.5|) the BGK approximation was used St(fx) = 
—vfi, where v is the transport collision frequency and introduced the new variable: the 
total energy along x-axis e x = + <p(x). There have been a number of studies, which 



explored the effects of the exact collision integral on collisionless phenomena [^TJ > 1 52| . These 
treatments use expansion in series of spherical functions in velocity spaces. The exact 
calculation are important only if the collision frequency is a strong function of the polloidal 
scattering angle. If the differential cross section does not depend on the polloidal scattering 
angle, the BGK approximation is correct exactly |j2]. For partially ionized plasma electron- 
neutral collisions are the most frequent scattering mechanism. At typical electron energies 
in the low-pressure discharges 1-5 eV |J, the differential cross section weakly depends on 
the polloidal scattering angle, and, therefore, the BGK approximation has a good accuracy. 

Equation ( |2.5| ) can be solved by a number of different methods. First, let us consider 
a direct solution. Alternative derivation using Fourier series is performed in Appendix C. 
After some straightforward algebra described in Appendix A, the symmetric part of the 
EVDF f ls = l/2(/ 1+ + is given by 



/u(v,x) = -mv y Vy f (x,e x )^, 



(2.6) 



where Vyf(x, e x ) = 1/2(V^+ + Vyl), Vy± are the oscillatory velocities of an electron with a 
given e x , ± signs denote v x > and v x < 0, respectively; 



m sinh $_ 



T 



cosh$ /*+ E y (x' 
cosh($+ - $) g_ 

dx 



vrf(„ e „ ^ _ e cosnv Jx ■ n y \x jcosh($+ - $')dr'+ , _ 

* •' • - ,i • • j„ i ,..,„,,, ,,> ^ r x E y (x') cosh &dr' ' 10 

(2.8) 

(2.9) 
(2.10) 

where x-(e x ), x + (e x ) are the left and right turning points, respectively, for the electron with 
energy e x [corresponding to zero velocity v x or e x = etp(x-)], r is the time of flight from 
the left turning point x_{s x ) to x, v± = Jv% + v%. The functions and $ depend on the 



V X [ X. Sr. 



rx 

$(x,e x ,v±)= {-iuj + v)dr, 

J X- 



electron speed via the collision frequency v{v). 

In the local limit the electron mean free path is large X » 5 and phase 7m($) >> 1. 
Therefore, cosh$ w sinh$ w l/2exp($). The main contribution in the both integrals in 
Eq.( [7.10|) are near the point x' — x and hence d$ = (—iu + v)dr 



yrf 
y 



Ey(x) 



m —iuo + v ' 



(2.11) 



as it should be in the local limit. 
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III. CALCULATION OF NONLOCAL CONDUCTIVITY 



Knowing the EVDF /i s , one can calculate the current density 

,3/2 



em J '' f , , 3 



(3.1) 



Substituting / ls from Eq. ( |2.6D into Eq.( |3.1| ) and making the transformation to the spherical 
coordinates in the velocity space dv x dv y dv z = v 2 dv sin 'dd'Odtfj in (cost? = v x /v; tan-?/> = v y /v z 
) Eq. (|3.1|) becomes 

j(x) = -ev^m 3 / 2 w (v 2 y V y rf ) ^^*>> ( 3 - 2 ) 
where the averaged over velocity direction factor (v y V y fj is 

(v 2 y V y rf ) = ^ [ J'* V y rf (x,e x ,v){sm<df{cos^ 2 d^d$. (3.3) 

Because V^does not depend on ip, the integration over ^-angle can be completed. Changing 
integral from i? to v x = v cos $ gives 



or 



(v 2 y V y rf ) = ^ £ ^ r/ (^, e a , ^)(^ 2 - 
V y f {x,e x ,v ± )- 



v V rf ) = 



--de. 



(3.4) 
(3.5) 



(p{x) 



Substituting Eq. (pT5|) into Eq.(|3.2|) and changing integration from v to e yields 



e 
2 Jo 



e - s a 



fix) ^£ X - V [ X ) 



V y rf de x 



dfoje) 
de 



de. 



(3.6) 



Further simplifications are possible if the collision frequency v is small [y << uo) or v does 
not depend on electron velocity. In this case V y f(x,e x ,v±) is the only function of (x,£ x ). 
Integrating Eq. (|3~6l ) in parts yields 



3\ x ) 



2 Jtp{x) 



V y rf (x,e x ) 
yje x - p(x) 



de. 



fo(e)de. 



(3.7) 



If Vy* is a constant Eq.(|3/7|) gives trivial result: j = enV % 



rf 



Introducing a new function T(e) 



T(e) = 



(3.8) 



and integrating Eq. ( |3.7|) in parts one more time gives 



- V;f(x,e)T(e) d£ 



2 -M*) ^Je - <f(x) 



(3.9) 
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For the Maxwellian EVDF /„, Eq. 
Substituting Eq.(|2.7|) into Eq. 



is equivalent to Lieberman's et al. result |34 . 
yields the nonlocal conductivity operator 



[ X G(x, x')E y {x')dx' + [ L G(x', x)E y (x')dx' 

JO Jx 



(3.10) 



where 



G(x,x') 



1 

2 



cosh <3> cosh ($ + - $') 



Vie) 



2m Jmax(ip,ip') 



sinh $_ 



^Je - ip(x)^/e - ip(x') 



de. 



(3.11) 



Note that G(x,x') has a logarithmic singularity at x = x' | 3^| , but because calculation of 
the electron current in Eq. (|3.10| ) requires additional integration, there is no singularity in 
the current. 

In the limit of large gap, where 5 < A << L, i?e($) >> 1 and cosh $ cosh($ + — 
$')/sinh$ + — > l/2[exp($ — $')+exp(— $-$')]. And the region of integration in Eq.( pTT0[) 
beyond the skin layer can be omitted. In the local limit, where A << 5 , Eq. fl3.10|) gives the 
standard local conductivity, see Eq. fl2.lip . 



IV. CALCULATION OF THE TRANSVERSE RF ELECTRIC FIELD PROFILE 



Maxwell's equations can be reduced to a single scalar equation for the transverse electric 



field [28 



d 2 E„ 



+ 



,2 E V 



\j(x) + IS(x) - 5 anti IS(x - L)\ 



(4.1) 



dx 2 " c 

where the electron current j is given by Eq. (|3.10|) , / is the current in the coil , 5 ant i = 0, 
if the there is no any coil with the current located at x = L, and 5 ant i,k = 1, if there is a 
coil with the current —I at x = L. The ID slab system of two currents flowing in opposite 
directions describes very well a cylindrical configuration, where a coil produces rf currents at 
both plasma boundaries x = -R and x = R, R = L/2 @, |£§. The Eq.(PD and Eq.flHOD 
can be solved numerically using a finite difference scheme. There is major difficulty in such 
approach. Straightforward computing the complex Green's function in Eq. ( |3.11| ) is slow and 
time consuming p7| . The better approach is to solve the integro-differential Eq.( |4.1| ) making 
use of a spectral method. In ID geometry the electric field can be represented as a sum of 
harmonic functions. The method is described in the Appendix D. 



V. AVERAGING OF KINETIC EQUATION FOR THE MAIN PART OF THE 
EVDF 

Kinetic equation for / averaged over the discharge period is 



dfo e df e 

Vx^— H E sc (x)- h —Re 

ox m dv„ 2m 



yj 



St(fo), 



st(f) = st v M) + suf) + st ee (f) + st inel (f), 



(5.1) 
(5.2) 
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st v M) 



st ee (f) 




(f - f)da, 




vdw 



(W el f) 



vD ee —f) + — 
vdw \ ee dw J vdw 



Stinel(fo) — ^ 



(w + e%) 



w 



vV ee f) 



(5.3) 
(5.4) 
(5.5) 

(5.6) 



where w = mv 2 /2 is the kinetic energy, St^(f) is the part of the elastic scattering collision 
integral with differential cross section da, which changes the electron momentum but does 
not alter the kinetic energy, St £ el (f) accounts for energy change in elastic collisions, St ee (f ) 
is the electron-electron collision integral, and Sti ne i(fo) is the sum over all inelastic collisions 
with the electron energy loss e* k and inelastic collision frequency v% (see details for ionization 
and wall losses in Hi 



Here, the coefficients D ee , V ee , V e i are given by ^4j, (49 



V r 



el 



2m 



-wv, 



VL = 2wv f 



Jq dw^/wf 



n 



4 r dww^f + w 3 / 2 r 



dwf 



4nA P 



n 



m 2 v 3 



(5.7) 



(5.8) 



(5.9) 
(5.10) 



where v ee is the Coulomb collision frequency and A ee is the Coulomb logarithm. Note that 



at large electron energies V ee = 2wv ee and D e 



2wT e v ee , where T e 



2/3 J™dww 3 / 2 f/n. 



If the electron mean free path is large compared with gap (A >> L), the first two terms on 
the left hand side are dominant. Therefore, fo is approximately a function of the e x only, not 
a function of both variables x, v x separately. Similarly, St^i(f) is the largest term from the 
remaining terms in the equation. Therefore, fo is approximately isotropic, and is a function 
of e only, so that St(f ) = |5T]| . This assumption was verified by comparison with particle- 
in-cell simulations in Ref. || for a capacitive coupled plasma, in Ref.|§, [[3^] for a inductive 
coupled plasma, and in Ref.|20| for a ECR discharge and in Ref.f)3| experimentally. 

To find fo, it is necessary to average over fast electron bouncing and over all velocity 
angles. First, let us average over fast electron bouncing. In order to do so, we integrate all 
terms of Eq.(|5.1|) over the full period of electron bouncing 



j> dtTerm(x, v x ) 



x + dx 



Term[x, v x > 0) + 



dx 



Term(x, v x < 0), 



where Term(x, v x ) is a term in Eq.( |5.1| ). Because the first two terms represent the full time 
derivative df/dt along trajectory, they disappear after integration, and Eq. ( |5.1| ) becomes 



dt Re 

2m 



E*(x)^- 



dtSt(f ) 



(5.11) 
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Second, we integrate Eq. (|5.11|) over all possible perpendicular velocities dv y dv z with a given 
total energy m[v 2 y + v 2 z ) < 2e P5| , |50[ . 
Total averaging is a triple integral 

Term(x, v) = — J J dv y dv z j> dtTerm(x,v), (5-12) 

where the factor 1/Att is introduced for the normalization purposes. Note, that integral 
Eq.( [5.12D is simply averaging over all phase space available for the electron with the total 
energy e. 

4~~ J I ^ Vy( ^ Vz _/ dtTerm(x, v) = — J dxd 3 v5 [e — w — (f(x)] . (5.13) 

If the Term(x, v) depends on velocity only only via speed v, like ,e.g., the inelastic collision 
integral Eq. fl5.6j ), then integration of Eq.( F3.13p in spherical coordinates gives 

rx+ 

Term(x,v)(e) = / dxv(x,e)Term[x,v(x,£)}, (5-14) 

v(x,e) = \j2[e — ip(x)]/m. (5.15) 

Thus, averaging of the terms, which are functions of the position and the kinetic energy 
only, reduces to the integrating over the entire available discharge volume weighted with the 
velocity, which is the standard procedure that also appears in the averaged kinetic equations 
for the local (collisional) case O . 



A. Calculation of the nonlocal energy diffusion coefficient 



The averaged energy diffusion term originates from the averaged left hand side of 
Eq. (|5.11 ), which gives 



eE y (x, t) dfi e 



m dv„ 



•>7T 



J dxd 3 v5 [e — w — ^p(x)] Re 



E*(x 



dfi 

dv. 



y J 



(5.16) 



Using chain rule for the integration in dv y and the fact that dS(e — w — if)/ dv y = mv y d8{e 
w — f)/de, Eq. (|5.16|) becomes 



eEyix^c^ = emd Re r _ E * {x)fi . (517) 

m dv y Sti de J y 

Substituting fx from Eq. (PDQ and integrating in the velocities v y and v z yields 



eE y (x,t) dfx d df 



m dv y 

where D £ is the energy diffusion coefficient 



de £ de 



D.- 



Am Jo 



- — Re / de v (e — e. 



dx 



—e;{x)v^{x,e x ). 



(5.18) 



(5.19) 



Equation (|5.19| ) is the general expression for the energy diffusion coefficient: in the limiting 
regime of the small mean free path (A << 5) it tends to the local limit, in the intermediate 
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pressure range (5 « A << L) Eq. (|5.19|) corresponds to the hybrid heating: electron motion 
in the skin layer is collisionless, but the randomization of the velocity kick acquired during 
a single pass through the skin layer occurs due to collisions in the plasma bulk, and in the 
opposite limit of the large mean free path (A > L) Eq.( p.l9|) describes collisionless heating 
(see Appendix B,C for details). If the collision frequency does not depend on the kinetic 
energy the direct substitution of Eq. ( |2TT| ) for gives 



De{e) 
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oo 

E 



de x \E yn (e x ) \ 



v 



where 



Ey n i^a 



1 

7T 



E y {6) cos (n0) dO 



Note that D e (e) in the last equation accounts for the bounce resonance ftb(e x )n 
the transit time resonance u = v/S, which corresponds to the maxima of Ey n (^£ x ) ■ 



(5.20) 

(5.21) 
uj and 



VI. SELF-CONSISTENT SYSTEM OF EQUATIONS 

In summary, the self-consistent system of equations for the kinetic description of low- 
pressure discharges accounting for nonlocal and collisionless electron dynamics contains: 
1. The averaged kinetic equation for fo reads 



d 



dfo 



d 



de de 



Vee + V el ] f = E 



VI. 



w + et 



-f (e + et)-v*J 



(6.1) 



where the upper bar denotes averaging according to Eq. (|5.14|) and D ee is given by Eq. (|5.9| ), 
by Eq. (P) V d by Eq. (P) , and D £ by Eq. d£H) or by Eq. (pD . 

2. The rf electric field is determined from the Maxwell Eq. Q4.1|) , where the electron 
current is given by Eq. (|3.10|) . 

3. The electrostatic potential is obtained using the quasineutrality condition 



rii(x) = / f {e)Je- <p(x)de, 

Jip(x) 



(6.2) 



where rii(x) is the ion density profile given by a set of number particles and ion momentum 
equations [p7[| . Eq.( |6~2|) is solved in the form of a differential equation [|| 



<fy = T scr,, d H n i{ X )] 

dx e dx 
where T^ cr (x) is the electron screening temperature 

1 f°° „ , , de 



(6.3) 



Me) 



2n(x) J v (x) ^ £ _ 

4. The power deposition can be computed as 

1 



(6.4) 



P x) 



-Re 



E*(x)j(x) 



(6.5) 
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Substituting Eq.( |3.6|) and changing the integration order, Eq.(|6.5|) becomes 



P = -v 7 ^ [°° DJe)^P-d£. (6.6) 

JO CLE 

Equation (|6.6|) can be used as a consistency check. 

This research was supported by the U.S. Department of Energy Office of Fusion Energy 
Sciences through a University Research Support Program. The author gratefully acknowl- 
edge helpful discussions with R. Davidson, Badri Ramamurthi and E. Startsev. 



VII. APPENDIXES 



A. Derivation of /i 

Direct integration of Eq. (|2.5|) yields 



f 1+ (x,v) = - ev y-^: 



fi-(x,v) = -ev y 



dfo 
de 



e -[*(*)-*(*0]£/ dr / + / 1+ ( x _) e -*(«) 



JX- 



(7.1) 



(7.2) 



where ± signs denote v x > and v x < 0, respectively, and for brevity we introduced 
Ey = Ey(x') and dr = dx/\v x \. The two constants fx-(xJ) are to be determined 

from the boundary condition at the turning points, and 



PX 

(x) = / (— iuj + u)dr. 

Jx- 



Due to continuity of the EVDF 



(7.3) 



(7.4) 



Substituting the boundary condition at the turning points Eqs.fl7.4j) into Eqs. (|7. 1[) and ( [7.2[ ) 
yields 

- r + e*-*'E y / rfr' + / 1 „(x_)e* + = f + e'^^Ey'dr' + / 1+ (x_) e -*+, 



or 



sinh <3>_ 



cosh($+ - &)Ey'dr' 



(7.5) 



where $ = $(2;), $' = $(2;'), and $+ = $(a; + ). /1 enters into the current calculation only 
as a sum fi + + /x_. Therefore, we compute /i s = l/2(/i + + from Eqs. (|7.1| ) and (|7.2j ) 

A- = - ew ^ cosh $ - £ sinh($ - $')^(^)^'} , 

substituting /i + (x_) from Eq.( |r.5|) gives 

fis = -mv y V y ' — , 



(7.6) 



(7.7) 
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where 



v m sinh $_ 



(7.8) 



cosh $ /*+ cosh($ + - ®)E y 'dT'- 
sinh$+ /*_ sinh($ - ®)E y 'dT ! 

Splitting the first term into two integrals = + J^ + , and accounting for the fact that 
cosh $ cosh($ + - $') - sinh $ + sinh($ - $') = cosh $' cosh($ + - $') (7.9) 



gives 



V rf = - — 

v m sinh $_ 



cosh $ /*+ £„' cosh($ + - &)di J + 
cosh($ + - $) £,/ cosh $'dr' 



(7.10) 



B. Diffusion coefficient in the energy space 

The equation for the energy diffusion coefficient 

„ e „ f e 1 , , r x + dx 
D £ = - — Re / de x (e - e x 
Am Jo 

has correct limits in local and nonlocal cases. 



E* y (x)V y rf(x,e 



(7.11) 



1. local limit A << L 



In the local limit of the small mean free path, substituting Vy* from Eq. fl2.lip into 
Eq.( [7.11|) gives 



D F 



-Re / den- (e — e n 



dx E*{x)E y {x) 



Am 2 ' Jo — - J x _ yj2{e x -tp)lm ~ iu + v 
Changing the order of the integration and accounting for the fact that 



de x 



e — E n 



m 2 Jo * ^2{e x -v)/m 



(7.12) 



D F = -Re I dx\E y \ 2 — - 

(ur + v l ) 



6 



which corresponds to the local limit pi, IT? 



(7.13) 



2. nonlocal limit 5 « A << L 



In the nonlocal limit, collisions during the electron motion in the skin layer are rare. 
Therefore, Vy± is simply velocity "kick" due the rf electric field. Recalling that Vyf(x,e) = 

{Vyl + Vy£)/2, the last factor in Eq.flmp can be written as 



1 e 

2 m j. 



dx 



E*y(x)V y rf (x,e x 



(7.14) 
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where § dr is an integral along the electron trajectory entering and leaving the skin layer, 
and AVyoc is the total velocity kick after a single path through the skin layer, and the angular 
brackets denote averaging over phases of the rf field. Eq.( |7.11| ) simplifies to 



- \ f £ de x (e-e x )(AV y 2 oo ), (7.ir,) 



In the limit of a uniform plasma Eq.(|7.15|) was proposed in Ref. [p6|, |32 



3. collisionless limit A >> L 



The energy diffusion coefficient Eq. (|7.11|) with substitution of V£* from Eq. (|7.10|) is de- 
termined by the following integral: 

Int = r x + _ E*(x)V;fdr = r x + _ E* y (x)dr 

cosh $ f£+ E y {x') cosh($+ - <S>')dr' + cosh($+ - $) /*_ E v {ff) cosh $'dr' . 
The term in brackets is cosh <3> + cosh $ £L cosh $'<ir / + 



sinh $_ 



cosh $ /*+ E y sinh $'dr' + sinh $ E y cosh $'dr' 



Therefore, Int 



cosh $ 



sinh $4 



i /*+ cosh Vdr' / x s + £*' cosh $'(2t' + Jnf 1, where 



/nil 



/*+ £*dr cosh $ /*+ sinh $'dr' + sinh $ /*_ ^ cosh $'dr' 



Integrating in 



2/ l Jx w 

parts gives Intl = f*+ sinh$dr [/*_ [ E y' E y + E y E y] cosh ®' dr ' 

In the nonlocal sinh $ = i sin co>r + vt cos cjt. Because the energy diffusion coefficient is 
determined by the real part of the integral and the real part of the phase is small (~ v ) 
Intl can be neglected. Therefore, 

1 /■£ r x + r x + 

D £ = — Re de x (e — e x ) coth $ + / E J cos lot dr \ El" cos lot" dr" , (7.16) 

where sinh*!'-^ = i sin loT + vT cos loT. Main contribution comes from the points, where 
loT = Tin and coth $ + = 7r S n =i ^(^^ — Tin). 



2 oo 



^ re 

Ef(e x ) = f XH£x) E y (x') cos ut' dr' 

Jx-(e x ) 



(7.17) 
(7.18) 



This corresponds to the pervious results [{47|. 



C. Alternative derivations in Fourier space. 

The direct calculation described in the previous sections are rather cumbersome. The 
alternative derivation can be done easier using Fourier series. 

It is convenient to introduce the variable angle of the bounce motion 

iisqn(v x ) f x dx 
0(x,e x )= * \ x) / r-j-^, 7.19 
T(e x ) Jx- \v x {e x )\ 
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where T is the half of the bounce period for the electron motion in the potential well <p(x), 
which is given by 

dx 



T{e x ) 



(7.20) 



\v x (e x )\ 

The bounce frequency for the electron in the potential well is fi&(e x ) = n/T(e x ). Utilizing 
angle variable, Eq. ( j2.5|) simplifies to become 



- iwfx + ^b^-\e x + v y eE y (6)^ = -ufi, 



(7.21) 



where ± signs denote v x > and v x < 0, respectively and flb{s x ) = 7r/T(e x ) is the bounce 
frequency in the potential well. 

We shall use Fourier series in variable 9: 



g(x,e x ) = g n exp(in9) 



(7.22) 



9n= 2^ 



g(9, e x ) exp (—initO) d6 



(7.23) 



Note that in the last integral, the region < 9 < tt corresponds to v x > , and the region 
— 7r < 9 < corresponds to v x < 0. Utilizing the Fourier series Eq. ( |7.23| ), the Vlasov 
equation becomes 

(inQ b - iuj + u)f ln 



TP *Ml 

ynVy de ' 



where 



1 
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E y {9) cos {pjff) d9 



(7.24) 
(7.25) 



Making use of Fourier series Eq. ( |7.22| ) , Eq. ( |7.24| ) gives 



fls [Xf £a 



dfo 



X. Er, 



e 
m 



iii,v v V v J (x,e x ) — 

E 



y y ^^ x> de' 
E yn cos[n9(x)} 



inQb — iu) + v 



(7.26) 



Eq.( 7.26 ) is the alternative form of Eq.( 7.1U| ). 

Substituting the function Vyf(x, e x ) from Eq. (|7.26| ) into Eq.( p.9|) gives the current density 



3\ x ) 



2 oo 

— y 



r(e) Eyn cos[n0(x)] 



V (x) J £ _ <p( x -) inQ b -iu + u 



de. 



(7.27) 



The averaged energy coefficient is given by Eq. ( |7. 1 1| ) . Substituting the function Vy*(x,e x ) 
from Eq.( 7.26 ) into Eq.( 7.1l| ) gives 



D £ = -Re / de x (e - e x ) / — E* y (x) £ 
4 Jo J x _ v x y „±^ r 



E yn cos[n9(x)} 
inQb — iuj + v' 



or 



2 oo 



D £ (e) = — y J o de x \E yn (e x )\ 2 (e - e x ) 



2 • ( 7 - 28 ) 

flbi^x) [£lb(£x)n -co] +v 2 

Note that Eq. (|7.28| ) is valid for any collision frequency, and Eq.( |7.17| ) is valid only for u « to. 
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D. Solving the Maxwell equations for the rf electric field using Fourier Series 

System that has an antenna at x = and a grounded electrode at x = L in the uniform 
plasma were studied in Ref.[[2lJ. The papers |23] and [[31]] considered a cylindrical-like system 
in the uniform plasma. Both papers used Fourier series to solve Maxwell's equations. Here, 
we shall generalize the procedure for a case of a nonuniform plasma. 

Similarly to the previous analysis, it is convenient to continue the rf electric field sym- 
metrically E y {x) = E y (—x) outside of the slab. Then, the electric field is given by Fourier 



series 



E y( x ) = 5I S sCOs(£; s x), 



(7.29) 



s=0 



where s is an integer, k s = (2s + 1)tt/(2L), for the case of the grounded electrode, and 
k s = (2s + l)vr/L, for the case of the cylindrical-like system. Substituting (|7.29|) into 
Eq.( f4.1|) and integrating with the weight cos(k s x)/L over the region [— L, L] yields 



—k. + 



Js + 



I [1 + 5, 



anti,k\ 



where 



2 r L 

js = y / j(x) cos(k s x)dx. 



LJo 



Substituting the equation for the current density Eq.( [7.27| ) gives 



ne 



m sQbT 



1=0 



io + iv 



sQ 



bT 



where QbT = vttc/L, and we introduced the generalized plasma dielectric function 



^7.30) 



(7.31) 



(7.32) 



7 gen 



10 + %V 



sQ 



bT 



2 sQ 



bT 



m L 



E 



Tie) 



TV 



o inVt b (e x ) - isVt hT ^ Vt b (e x ) 



G s ,n( £ )Gi t n(e)de, 



(7.33) 

where the coefficients G^ n (e) are the Fourier transform of cos(kix) in the bounce motion of 
the electron in potential well: 



Gi >n (s) - Tj; 



cos[k[x(r)] cos 



/ 7mT\ 



dr 



V T J 

In the limit of a uniform plasma t = x/v x , r = L/v x which gives 



Gi, n {e) = - 



J cos(kix) cos y~~r^j d x 



For a cylindrical-like system coefficients G^ n (e) are particular simple 

1 



G,Je) 



-Al\,\n\, 



(7.34) 



(7.35) 



(7.36) 
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and the generalized plasma dielectric function is 



z 9 s T' 



\0 = 6 Sjl Z(£), 
where Z(£) is the "standard" plasma dielectric function 
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-1/2 



oc^expj 

-oo t - 



(7.37) 



(7.38) 



Eq.( [7.32|) is identical to the results of Ref.|23j for cylindrical-like configuration uniform 
plasma with a Maxwellian EVDF. 

Coefficients Gi )Tl (e) can be effectively computed using the fast Fourier transform |55 



The off-diagonal coefficients are generally very small, that is why using this spectral method 
makes computing much faster than the straight forward finite difference method used in 
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